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Abstract 

The distribution of inclusion-rich domains in membranes with active two-state inclusions is 
studied by simulations. Our study shows that typical size of inclusion-rich domains (L) can be 
controlled by inclusion activities in several ways. When there is effective attraction between state- 
1 inclusions, we find: (i) Small domains with only several inclusions are observed for inclusions 
with time scales lO^'^s) and interaction energy [~ ©(keT)] comparable to motor proteins, (ii) 
L scales as 1/3 power of the lifetime of state-1 for a wide range of parameters, (iii) L shows a 
switch-like dependence on state-2 lifetime fc^^^. That is, L depends weakly on ki2 when ki2 < 
but increases rapidly with ki2 when ki2 > A;^2) the crossover /c^g occurs when the diffusion length of 
a typical state-2 inclusion within its lifetime is comparable to L. (iv) Inclusion-curvature coupling 
provides another length scale that competes with the effects of transition rates. 

PACS numbers: 87.16.Dg, 05.40.-a, 05.70.Np 
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INTRODUCTION 



Biological membranes are composed of lipids, proteins, carbohydrates and other materi- 
als Many experimental as well as theoretical studies have shown that the distribution 
of these molecules in a membrane is nonuniform. Instead, dynamical patch-like structures 
with typical size ranging from tens of nanometers to about one micrometer exist in a mem- 
brane [2]. The mechanism for the formation of these heterogeneous structures has been a 
major research topic in membrane biophysics in recent years. Possible mechanisms include 
(i) continuous membrane recycling may result in nonequilibrium domains in biological mem- 
branes [sl, and (ii) the coupling between lipid density and local membrane curvature provides 
a length scale for membrane domains ^. The most important conclusion from these stud- 
ies is that finite-size domains in membranes have to be sustained either by nonequilibrium 
processes such as recycling or by effective long-range interactions due to membrane elasticity. 

Another important trend in recent biomembrane research is to treat biomembranes as 
nonequilibrium systems j^, 1^, Q, y, [ol, 1^, 
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which contain active inclusions that are 



driven by external stimuli. Earlier work on these active membranes studied mainly the effect 
of nonequilibrium activities on their fluctuation spectrums Q, 



to activities were also discussed in 
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12[ |. the instabilities due 



13| . In these models, the effect of conformation change 
of the inclusions are neglected. On the other hand, it has been shown by Lacoste and Lau 6| 
that active inclusions with an internal time scale have different fluctuation power spectrum 



from simpler mode^ 
were predicted in 



s which neglects this effect. Furthermore, finite-wavelength instabilities 
5, 



, Is, Q] for membranes with two-state active inclusions, and it has 
been pointed out in [5] that these instabilities are similar to the dynamic nm-scale domains 
in biomembranes. Indeed, recent experiments have observed enhanced spatial aggregation 
of bacteriorhodopsin (BR) due to light-induced activities in model membranes IJ] and 
clustering of light-harvesting antenna domains in photosynthetic membranes under low- 



light conditions pj|. These experimental observations suggest that size of inclusion-rich 
domains in a membrane can be tuned by changing the strength of external stimuli of the 
active inclusions. Therefore in this article we study the physics of activity-controlled finite 
size domains in membranes in a simple model. 

To focus on the basic interactions responsible for inclusion-activity induced finite size 
membrane domains, we consider the simplest possible model, i.e., a fluid membrane com- 
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posed of one type of lipid and two-state active inclusions. An inclusion is an active unit which 
corresponds to a molecular complex that contains proteins, lipids, and other molecules. As 
shown schematically in Fig. 1, inclusions in different internal states have different confor- 
mations, therefore the interactions between the inclusions and lipids depend on the internal 
states of the inclusions. There are also couplings between the inclusion densities and mem- 
brane curvature due to the up-down shape asymmetry of the inclusions. Lattice Monte 
Carlo simulations are applied to study the steady-state distribution of the inclusions when 
there is attractive interaction between state- 1 inclusions. The effect of inclusion-curvature 
couplings and inclusion transition rates kajj (probability that a state-/? inclusion transforms 
to state a per unit time) on the typical size of inclusion-rich domains is discussed. The main 
results of our simulations are: (i) Small domains with only several inclusions are observed 
for inclusions with time scales (~ 10~^ s) and interaction energy [ ~ (9(kBT)] comparable to 
motor proteins, i.e., our simple model is able to produce the kind of small inclusion clusters 
observed in [l41] within typical parameter range, (ii) Typical size of inclusion-rich domain 
(L) scales as L ~ k2i^^^'^. This often observed behavior 0, Q, [sl provides a mechanism 
to control continuously the size of inclusion-rich domains by inclusion activities, (iii) A 
switch-like response of typical size of inclusion-rich domains to the external stimuli is found 
at fixed k2i: There exists a crossover transition rate /cjfg such that L depends weakly on the 
stimuli (ku) when ki2 ^ but becomes very sensitive to it (i.e., L increases rapidly as fci2 
increases) when ki2 ^ kl2- The crossover A;^2 occurs when the diffusion length of a state-2 
inclusion within its lifetime is about the same as L. (iv) L decreases when the coupling 
between state-1 inclusion and membrane curvature increases. 



MODEL 



The system is described by a lattice model with membrane height h{i,j), inclusion den- 
sities (paii.i) (a = 1 or 2), and lipid density (f)Q{i,j) defined on a two-dimensional square 
lattice. The Hamiltonian of the system H = Hm + Hi + includes the elastic energy of 
the membrane, the interaction energy between the inclusions and lipid molecules, and the 
inclusion-membrane coupling. To lowest order 

= [«:(ViM^,j))' + 7(VxM^,j))'] , (1) 
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where k is the bending rigidity, 7 is the surface tension of the membrane, a is the lattice 
constant, V]_h{i,j) = [h{i + 1, j) + h{i — 1, j) + h{i,j + 1) + h{i,j — 1) — 4/i(z, j)]/a^, and 
V±h{i,j) = i[h{i + 1, j) - h{i - 1, j)]/2a + ][h{i, j + 1) - h{i,j - l)]/2a. Hi includes the 
short-range interactions between the inclusions and the lipids, the simplest form is nearest- 
neighbor interaction 

2 

Hi= ^ Jmn(pmii,j)<Pnik,l), (2) 

{{i,j){k,l)) ni,n=0 

where Ylm=o^rn{'>',j) = 1 for all '^((ij)(ki)) means sum over all nearest-neighboring 

pairs. The simplest form for the coupling between the inclusion density and local membrane 
curvature is 

2 

Hc = ^^ KCm(j)m{hj)[Vlh{i,j)], (3) 

where cq is the spontaneous curvature of the lipids, Cq, (a = 1, 2) is the coupling between 
state-a inclusions and local membrane curvature. The inclusion-curvature coupling constant 
Ca ~ ASa/Sa/o, has dimension of inverse length, where S„ is the average cross sectional 
area of a type-a inclusion, AEq is the difference between head area and Sq,, and /„ is the 
thickness (in the z direction) of the inclusion [l^. 

The simulations are performed on a Lm x Lm lattice with periodic boundary condition. 
The lattice constant a ~ 5 nm is chosen to be the smallest length scale for the continuum 
elasticity theory of membranes to be valid [17]. In the simulations each lattice site is either 
occupied by a state-1 inclusion, a state-2 inclusion, or lipids. That is, (j)m{i,j) = or 
1 (m = 0, 1, or 2), and Ylm=o ^rn{h j) = 1 for all Since our goal is to study 

the effect of kaf^ and Cm on the distribution of the inclusion-rich domains, in the present 
study the active forces exerted by the inclusions and the viscous flow of the solvent are, 
similar to another active membrane simulation 18|, neglected. Each Monte Carlo step 



is composed of three independent trial moves: the motion of the membrane, the in-plane 
motion of the inclusions, and the conformation change of the inclusions. The dynamics 
of h{i,j) is simulated by updating the membrane height of a randomly chosen lattice site 
with a displacement between —Ah and Ah with Metropolis algorithm, where Ah = 0.4 
nm. The in-plane motion of the inclusions is simulated by Kawasaki exchange dynamics. 
Time interval of one Monte Carlo step can be derived from the two-dimensional diffusion 



law (r^) = ADt. Since typical diffusion constant for membrane proteins is on the order 
of Ijjm^/s [3], and a free inclusion moves a distance a in one Monte Carlo step, a single 
Monte Carlo step corresponds to time interval At ~ lO^^s. The inclusion conformation 
changes are simulated by assigning each state-/? (state-a) inclusion a probability kapAt 
{kpoL^t) to change its state to state a (state /?) in each Monte Carlo step. We choose kaps 
to be constants in the simulations, i.e., we assume that the coupling between conformation 
change of the inclusions and local membrane curvature and composition is small. Therefore 
the conformation change does not obey detailed balance, the steady state of the system is 
not in thermal equilibrium 20|. In addition, in this article all the simulations are performed 
with kapAt < 0{1) such that the lifetime of each state is sufficiently long that inclusion 
diffusion can be taken into account properly. The initial condition of the simulations is 
randomly distributed A^inc = 0inc x Lm x Lm inclusions with A^inc/2 of them in state 1 and 
of them in state 2 in a flat membrane. The probability that an arbitrary inclusion 
is within an inclusion cluster of M sites, P{M), is analyzed when the systems have reached 



steady states. Typical size of inclusion-rich domain is taken as L = i/Mmax, where Mmax is 
the value of M which maximizes P{M). 

The free parameters in our model include c^, Jmn {m = 0, 1, 2), and ki2, ^21- Jmn arises 
from short range noncovalent interactions between the molecules including the effect of 
hydrophobic length mismatch between different molecules or inclusions of different internal 
states. Since we are interested in the case where finite-size inclusion-rich domains exist 
in the steady state, Jmn are chosen such that state-1 inclusions attract each other, and 
state-2 inclusions have very weak or no attraction with other inclusions and lipids. To 
prevent aggregation of state-2 inclusions due to membrane curvature-induced attraction, 
C2 has to be close to Cq. For state-1 inclusions, there is no such restriction on the choice 
of Ci. The simplest choice of parameters that satisfies the above criteria is Cq = C2 = 0, 
ci > 0, Joo = Ju = J22 = J02 = J, and Jqi = J12 = J + AJ with AJ > 0, and this 
choice is applied to our simulations. The bending elastic modulus and surface tension of 



the membrane are chosen to be typical values, k = 5 x 10^^° N-m 2l|, and 7 = 24 x 10~^ 
N-m ^22]. A J ~ (9(kBT) because the interactions between inclusions and lipids are non- 
covalent, typically Ci ~ ASi/Si/i ^ 0{1 nm~^). For convenience we define a dimensionless 



constant G = C\a^ ksT j k to describe the inclusion-curvature coupling, its typical value is 
G<0{\). 
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The transition rates kaf3 for the inclusions depend on the mechanism of conformation 
transitions of the inclusions under consideration, (i) For inclusions have an intrinsic time 
scale (like ion pumps), an inclusion stays at state 1 (2) in the absence of stimuli, and 
changes its conformation to state 2 (1) when it binds to energy source such as ATP or 
specific ions, the lifetime of state 2 (1) is an intrinsic property of the inclusion, it is on the 
order of 10^^ — 10~^ s, therefore ku {k2i) ~ 10^ — 10^ s^^. (ii) For inclusions that change 
conformation by binding to or dissociating with small ligands, the transition rates ki2, ^21 
can exceed 10^ s~^, depending on the concentration of ligands in the solvent. 



DISCUSSION 



Figure 2 shows domain size distribution P{M) for 0inc = 12.5%, AJ = l.bksT, 
ki2 = 10^ s~^, k2i = /i;i2/32, and G = 0, 1, and 2. The peak at M = 1 comes from the isolated 
inclusions in the inclusion-poor domain. The peak at greater M provides the characteristic 
size of inclusion-rich domains. As G increases, the characteristic size of inclusion-rich do- 
mains decreases and the peak of P{M) becomes more significant due to inclusion-membrane 
coupling. This effect can be better visualized by the snap shots from the steady states. As 
shown in Fig. 3, when G ^ the location of inclusion-rich domains has strong correlation 
with the regions with high membrane curvature. Since the system has lower free energy 
when the inclusions are in the regions with higher curvature, thus the membrane forms 
many mountain-like regions with inclusions in them. This effect makes the width of P{M) 
decrease as G increases and the position of the peak shifts to smaller M. 

Figure 4 shows the relation between L and k2i for 0inc = 12.5%, AJ = l.dksT, and 
ki2 = 10~^ s~^ on a 128 x 128 lattice. The choice of ki2 is made to study inclusions with an 
intrinsic time scale ~ 10 ms, or inclusions in a solution with high density ligands that induces 
2 — s> 1 transitions. As long as AJ is order unity and 0inc is not too small for the simulation 
system size, simulations with other choices of A J and (p^^^ give results that are qualitatively 
the same as Fig. 4. First of all, typical size of inclusion-rich domains agrees with L ~ ^21"^^^ 



pretty well for a wide range of k2i, and the agreement is better when G is smaller 



23j. This 



inclusion activities is a two- 



is because the growth of inclusion-rich domains in the absence o 
dimensional phase separation dynamics, i.e., L ~ t^/^ |2j, y, 
saturates due to state-1 to state-2 transitions on time scale ~ k2i~^, thus the typical length 



26| . This growth eventually 
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scales of inclusion clusters in the steady state for G ^ 0{1) agrees with L ~ ^21"^^^ [s^. 
This relation is also observed in simulations with other 0inc, Jmn-, and ki2 (data not shown). 
Finite size domains with 1/3 scaling law due to similar mechanisms have been observed in 
domains in membranes induced by continuous recycling processes [3| , and membranes made 
of binary reactive lipids {3, Q] ■ Another interesting result is that small inclusion-rich clusters 
with Mjnax ~ 10 are observed when both A;2iAt > 10~^ (i.e, /c2i ~ 10^ s"-*^) and ki2^t < 10~^ 
(i.e., ki2 ^ 10^ s^^) are satisfied. Therefore in systems with interaction energy and active 
transition rates comparable to the characteristic time scale of a motor protein, the lifetime 
of state- 1 inclusions could be too short for a large cluster to be formed. In enhanced 
clustering of light-sensitive bacteriorhodopsin (BR) is observed in the presence of light, but 
there are only a few BR in a typical cluster. A possible explanation for this observation 
is that in the absence of light, BRs have very weak interaction with each other (like our 
simulation with ki2 = and all inclusions are in state-2); while in the presence of light, BRs 
in an intermediate state (analogous to state- 1 inclusions in our simulation) have attractive 
interaction with each other, so BRs begin to form clusters. However, the lifetime of this 
intermediate (for BR, it is on the order of 1 ms) is too short for large clusters to be formed, 
thus only very small BR clusters are observed. 

Figure 5 shows the relation between L and ki2 for 0inc = 12.5%, A J = l.bksT, and 
k2i = 10 s~^ on a 128 x 128 lattice. Here L shows a switch-like dependence on ki2: L depends 
weakly on ku when ki2 is less than some crossover rate fcjj'g but increases rapidly with ku 



when ki2 ^ ^12- This can be understood by comparing L2 ~ y [^) ki2~^, the diffusion 
length of a state-2 inclusion within its lifetime, to L. When L2 ^ L, state-2 inclusions in 
inclusion-rich domains can escape to inclusion-poor domains within their lifetime, thus L 
depends weakly on ki2. On the other hand, as illustrated in the inset concentric circles of 
Fig. 5, when L2 ^ L, the chance that a state-2 inclusion in an inclusion-rich domain cannot 
escape to inclusion-poor domains increases rapidly as AL = L — L2 increases. In this case 



L increases rapidly with ki2. Thus fcjfg is determined by the condition y (-^j A;*2^^ ~ L. 
As can be seen from Fig. 5, for given G, the crossover occurs when ki2At ~ (L/a)~^, and 
the switch-like behavior is more significant for systems with smaller G. For simulation 
with k2i > O{10 s~^), smaller L causes fc*2 to be larger, and the switch-like behavior is 
less significant because simulations are performed with kuAt < 0{1). In the limit of very 
large ki2, one expects that almost no state-2 inclusion can escape from the domain within 
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its lifetime, and typical size of inclusion-rich domains should saturate to its equilibrium 
value either determined by system size or (when G is large) by the length scale selected 
by inclusion-membrane coupling. However, this limit happens when ki2At > 0{1), and 
therefore cannot be accurately simulated in the present work. 

The inclusion-curvature coupling G also has important effects on L. First, both Fig. 4 and 
Fig. 5 show that L decreases as G increases. That is, the coupling between inclusion density 
and membrane curvature suppresses the formation of large inclusion clusters because the 
inclusions tend to locate in regions with greater membrane curvature. As a result, instead 
of forming large inclusion-rich domains, the membrane prefers to form many mountain-like 
regions with inclusions aggregating on them when G ^ l'^. Second, L ~ k2i~^^^ relation 
in Fig. 4 does not describe large G well as small G cases because besides the length 

scale set by inclusion transitions, inclusion-curvature coupling provides another length scale 
for inclusion- rich domains. This is also true in Fig. 5, where the rapid increase of L at 
ki2 > ~ {L/a)~^/At is not significant at large G. Thus, in general inclusion-curvature 
couplings competes with the effects of kap- 

In summary, our study shows that there are several ways for inclusion activities to control 
the typical size of nonequilibrium membrane domains. In particular, the power-law relation 
L ~ k2i~^^^ provides a way to continuously tune the size of domains, the switch-like de- 
pendence of L on ki2 provides a mechanism for sudden change of domain size dependence 
on the inclusion activities. Although our model has neglected hydrodynamics of the solvent 
and active forces exerted by the inclusions, the mechanism proposed in the current study is 
rather general and should exist in models which take these effects into account. A detailed 
analysis that includes the effects of active forces and full hydrodynamics will be presented in 



a future work 



27| . The small inclusion clusters in the large /c2i regime observed in our sim- 



ulations are similar to the experimental observations in We expect more experiments 
on these systems and other membranes can be performed to study the mechanisms revealed 
in our model. For example, the different hydrophobic lengths of different internal states of 
the rhodopsin 28| can be used to induce effective interaction between nearby rhodopsin, 
therefore a rhodopsin-containing membrane is also a good candidate for experimental study. 

This work is supported by National Science Council of the Republic of China under Grant 
No. NSC-93-2112-M-008-020. Part of this work is done during HYC's visit to the National 
Center for Theoretical Sciences (NCTS), Hsinchu, Taiwan. 
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ure captions 

» FIG. 1: Schematics of a membrane with two-state inclusions. 

• FIG. 2: Domain size distribution P{M) for ^i^c = 12.5%, A J = IMbT, ku = 
10^ s~^, k2i — A;i2/32, and G = 0, 1, and 2. The peak at M = 1 comes from the 
isolated inclusions in the inclusion-poor domain. The peak at greater M provides the 
characteristic size of inclusion-rich domains. As G increases, the characteristic size of 
inclusion-rich domains decreases and the peak of P{M) becomes more significant due 
to inclusion-membrane coupling. 

» FIG. 3: Snapshots for steady state inclusion distribution (left, dark regions are oc- 
cupied by inclusions) and membrane height (right) for (f)^^^ — 12.5%, AJ = 1.5A;bT, 
fc^2 = 10^ s-\ and k2i = ki2/?>2. (a). G = 0, (b). G = 2. In order to faithfully 
present the morphology of the membrane, the unit length for h and the unit length 
in the xy plane are both chosen to be a. For nonzero G membrane curvature close 
to the inclusion-rich domains increases and the typical size of inclusion-rich domain 
decreases. 

I FIG. 4: Typical size of inclusion-rich domains in the steady state for 0inc = 12.5%, 
AJ = LbkeT, ku^t = 10"=^. G = (circles), G = 0.5 (squares), G = 1.0 (diamonds), 
and G — 2.0 (triangles). The dashed hne has slope —1/3. 

» FIG. 5: Typical size of inclusion-rich domains in the steady state for ^inc = 12.5%, 
AJ = 1.5A;Br, k2iAt = 10"! G = (circles), G = 0.5 (squares), G = 1.0 (diamonds), 
and G — 2.0 (triangles). Data with greatest value of ki2A.t are obtained from simu- 
lations on a 256 x 256 lattice. The inset is a schematic of a domain with radius L. 

When ki2 > /i;*2, state-2 inclusions in the inner circle (radius AL) have small chance 
to leave the domain within its lifetime. AL = L — L2 increases as ki2 increases. 
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FIG. 2: Domain size distribution P{M) for (^i„c = 12.5%, AJ = 1.5/csT, ki2 = 10^ s-\ k2i = 
/ei2/32, and G = 0, 1, and 2. The peak at M = 1 comes from the isolated inclusions in the 
inclusion-poor domain. The peak at greater M provides the characteristic size of inclusion-rich 
domains. As G increases, the characteristic size of inclusion-rich domains decreases and the peak 
of P{M) becomes more significant due to inclusion-membrane coupling. 
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FIG. 3: Snapshots for steady state inclusion distribution (left, dark regions are occupied by in- 
clusions) and membrane height (right) for (/)inc = 12.5%, A J = l.^ksT, ki2 = 10^ s^^, and 
^21 = A;i2/32. (a). G = 0, (b). G = 2. In order to faithfully present the morphology of the 
membrane, the unit length for h and the unit length in the xy plane are both chosen to be a. For 
nonzero G membrane curvature close to the inclusion-rich domains increases and the typical size 
of inclusion-rich domain decreases. 
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FIG. 4: Typical size of inclusion-rich domains in the steady state for ^inc = 12.5%, AJ = LSfesT, 
ki2^t = 10~^. G = (circles), G = 0.5 (squares), G = 1.0 (diamonds), and G = 2.0 (triangles). 
The dashed line has slope —1/3. 
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FIG. 5: Typical size of inclusion-rich domains in the steady state for ^inc = 12.5%, AJ = l.SfcsT, 
fc2iAt = 10^^. G = (circles), G = 0.5 (squares), G = 1.0 (diamonds), and G = 2.0 (triangles). 
Data with greatest value of k\2^t are obtained from simulations on a 256 x 256 lattice. The inset 
is a schematic of a domain with radius L. When k\2 > ^12 > state-2 inclusions in the inner circle 
(radius AL) have small chance to leave the domain within its lifetime. AL = L — L2 increases as 
ki2 increases. 
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